Method for improving the spatial resolution of a compton camera

ABSTRACT

A Compton camera has two detectors which each sense the location, energy and time of photon collisions. Compton events are detected in which a photon emanating from the subject collides with one detector at energy (E 1 ) and then the second detector at (E 2 ), both energies E 1  and E 2  are measured and used in conjunction with an initial photon energy E 0  to increase the resolution of the camera.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a 371 of PCT/US98/09921 filed May 15, 1998 whichclaims benefit of Provisional No. 60/046,774 filed May 16, 1997.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

This invention was made with United States Government support UnderGrant No. CA-32846 awarded by the National Cancer Institute. TheGovernment has certain rights in the invention.

BACKGROUND OF THE INVENTION

The field of the invention is x-ray and gamma-ray cameras based on theCompton-scatter principle.

Scintillation cameras are widely used as diagnostic tools for analyzingthe distribution of a radiation-emitting substance in an object understudy, such as for the medical diagnosis of a human body organ. Atypical scintillation camera of the Anger-type is described in U.S. Pat.No. 3,011,057.

The Anger-type scintillation camera can take a “picture” of thedistribution of radioactivity throughout an object under investigation,such as an organ of the human body which has taken up a diagnosticquantity of a radioactive isotope. A scintillation camera of theAnger-type produces a picture of the radioactivity distribution bydetecting individual gamma rays emitted from the distributedradioactivity in the object that pass through a collimator to produce ascintillation in a thin planar scintillation crystal. The scintillationis detected by a bank of individual photomultiplier tubes which viewoverlapping areas of the crystal Appropriate electronic circuitstranslate the outputs of the individual photomultiplier tubes into X andY positional coordinate signals and a Z signal which indicates generallythe energy of the scintillation event and is used to determine whetherthe event falls within a preselected energy window. A picture of theradioactivity distribution in the object may be obtained by coupling theX and Y signals which fall within the preselected energy window to adisplay, such as a cathode ray oscilloscope which displays theindividual scintillation events as spots positioned in accordance withthe coordinate signals. The detection circuitry typically provides forintegrating a large number of spots onto photographic film.

The “resolution” of a scintillation camera refers to the ability of thecamera faithfully to reproduce the spatial distribution of theradioactivity which is within the field of view of the device. Theoverall intrinsic resolution of the Anger camera detector is generallydependent on the ability of the detector to signal accurately theposition coordinates of each scintillation event. In the Anger-typecamera there is no way to determine the direction of the gamma ray whichproduced the scintillation event merely by detecting the location of theevent. Instead, a collimator is placed in front of the detector torestrict to a narrow angle the direction of the gamma rays reaching thedetector.

As a result of the need for a collimator to restrict the angle of theincident gamma rays, the Anger camera has a relatively low sensitivity,which is defined as that fraction of the gamma rays emanating from thesource which actually reach the detector to produce an event thatcontributes to the image. If the collimator is less restrictive to admitmore gamma rays and increase sensitivity, the spatial resolution of thecamera is reduced.

Another type of camera relies on the Compton effect to determine theangle of incidence of a gamma ray. This camera does not require acollimator in order to determine gamma ray angle of incidence and a veryefficient camera is, therefore, possible. The operating principle of aCompton camera is described by Martin JB, Knoll GF, Wehe DK, Dogan N,Jordanov V, Petrick N: A Ring Compton-scatter Camera For Imaging MediumEnergy Gamma Rays, IEEE Trans. Nucl. Sci., 1993; 4-:972-978; and SinghM: Electronically Collimated Gamma Camera For Single Photon EmissionComputed Tomography, Part I: Theoretical Considerations and DesignCriteria, Med. Phys., 1983; 10(4) :421-427.

FIG. 1 shows a scatter aperture for a Compton camera which is composedof two detectors D1 and D2. Incident X- or γ-ray photons interact in thefirst detector D1 via Compton-scattering, and the scattered photon isdetected in time-coincidence in the second detector D2. Using themeasurements of the energy deposited in the two detectors, and locationof D1 and D2, the direction of the incident photon can be resolved towithin a conical ambiguity. Alternatively, if the energy of the incidentphoton is known, a measurement from only one of the two detectors isneeded. By collecting a large number of these interactions andsubsequent data processing, the conical ambiguity can be reduced and thesource intensity of incident photons recovered. The primary attractionof Compton-scatter apertures for single-photon imaging is that theypotentially offer reduced counting noise or improved spatial resolution.The potential for this improvement exists because the solid-angleefficiency of a scatter-camera can be two orders of magnitude higherthan a conventional mechanically-collimated system for equal spatialresolution.

Despite the large gains in raw counting efficiency, some advantage islost in the process of recovering the source distribution by removingthe intrinsic ambiguity. Further reductions in performance—especiallyfor imaging photons having energies below 200 keV—can arise fromdetectors having modest energy resolution. In a Compton-scatter camera,poor energy resolution translates directly to poor spatial resolution(and large noise amplification if this blurring is unfolded).

In mechanical collimated systems, efficiency must always be traded offagainst spatial resolution. High resolution necessarily means poorsensitivity. This is not the case for a Compton camera. Any method ofreducing the uncertainty in the direction of the scattering angle thatdoes not reduce efficiency will improve resolution. Using detectorshaving better energy resolution is one way to accomplish this.“Conventional” Compton cameras employ one detector having good energyresolution and another with only modest resolution. Measurements fromboth detectors are typically used to determine the incident energy E₀.Once the incident energy has been estimated, the measurement from thedetector having the best energy resolution (E₁ or E₂) is used toestimate the scattering angle.

BRIEF SUMMARY OF THE INVENTION

The present invention is an improvement for a Compton camera whichincreases its spatial resolution by reducing the uncertainty of thescattering angle. More particularly, the present invention is a camerawhich includes a first detector (D1) for detecting the time, locationand energy (E1) of photon collisions; a second detector (D2) fordetecting the time, location and energy (E2) of collisions with photonsemanating from the first detector (D1); a Compton event detector; and anangle detector for determining the angle of Compton scattering based onthe energies E₁, E₂ and an estimation of energy E₀ of the photonsemanating from the subject. The location and angle information for eachCompton event is used by an image reconstruction means to produce animage.

A general object of the invention is to improve the spatial resolutionof a Compton camera. The present invention reduces the angularuncertainty of the incident photon by using two detectors having highenergy resolution in conjunction with knowledge of the energy spectrumof the incident photons. Where the incident energy spectrum is knownexactly, the present invention can reduce the uncertainty in thescattering angle by up to a factor of {square root over (2)} over thebest single measurement. The maximum reduction occurs when two detectorsD1 and D2 having equal energy resolutions are employed. If detectorshaving unequal energy resolutions are used the uncertainty σ_(E) isgiven by $\begin{matrix}{\sigma_{E\quad} = \sqrt{\frac{\sigma_{1}^{2}\sigma_{2}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}}} & (1)\end{matrix}$

where σ₁ ² and σ₂ ² are the variances in the energy uncertainties forthe first and second detectors, respectively. While {square root over(2)} may seem like a small gain in spatial resolution, the pixel-wisevariance in reconstructed images, as a function of spatial resolution,can be quite non-linear. Small resolution changes can, therefore, resultin relatively larger variance changes in the intensity estimate at eachimage pixel.

In one embodiment the camera includes a third detector (D3) fordetecting the time, location and energy (E3) of collisions between saidthird detector (D3) and photons emanating from the second detector (D2);and the angle determining means is also responsive to the energy (E3) tocalculate the angle of Compton scattering.

In one embodiment the angle determining means calculates the angle φ ofCompton scattering according to the expression:$\left( {\hat{\phi},{\hat{E}}_{0},{\hat{X}}_{1}} \right) = {\underset{({\phi,E_{0},X_{i}})}{argmax}{f\left( {\phi,E_{0},{X_{1}\left\{ Y_{i} \right\}},\left\{ t_{i} \right\},\left\{ X_{i} \right\}} \right)}}$

where X₁ is the point of impact on the first camera andf(φ,E₀,X_(i)|{Y_(i)},{t_(i)},{X_(i)}) is given by equation:$C \times {\int{\cdots {\int{\left\lbrack {\prod\limits_{i = 1}^{N}\quad {{f\left( {Y_{i}E_{j}} \right)}{f\left( {t_{i}\tau_{i}} \right)}}} \right\rbrack {f\left( {\left\{ E_{i} \right\},\left\{ \tau_{i} \right\},{\left\{ X_{i} \right\} E_{0}},\phi} \right)}{E_{i}}{\tau_{i}}}}}}$

In another embodiment the angle determining means chooses an incidentenergy estimate E₀ as Y₁+Y₂ where Y₁ and Y₂ are measured energiescorresponding to E₁ and E₂ and the determining means calculates theangle of the Compton scattering by first estimating energy E₁ usingenergies Y₁ and Y₂ and estimate E₀ and then using the E₁ estimate todetermine the scattering angle φ.

The determining means may preferably estimate energy E₁ by solving thefollowing equation:

Ê₁=(σ_(D1) ²+σ_(D2) ²)⁻¹[σ_(D2) ²Y₁+σ_(D1) ²(Ê₀−Y₂)].

Wherein Doppler broadening (explained below) is negligible thedetermining means preferably estimates angle φ by solving the followingequation:${\hat{\phi} + {\cos^{- 1}\left\lbrack {\frac{\alpha}{{\hat{E}}_{0}} - \frac{\alpha}{{\hat{E}}_{0} - {\hat{E}}_{1}} + 1} \right\rbrack}}\quad$${\text{where:}\quad \alpha} = \frac{E_{0}}{m_{0}c^{2}}$

and where m₀c² is the rest mass of an election.

Where Doppler broadening is significant the determining means calculatesthe angle of the Compton scattering by solving the following equation:$\hat{\phi} = {\underset{\phi}{argmax}{{f\left( {{\hat{E}}_{0},{{\hat{E}}_{1}\phi}} \right)}.}}$

The invention also includes a method for use with the above referencedcamera wherein the method comprises the steps of detecting the time,location and energy (E1) of collisions between said first detector (D1)and photons emanating from the subject, detecting the time, location andenergy (E2) of collisions between said second detector (D2) and photonsemanating from the first detector (D1), identify pairs of collisions inthe respective first and second detectors that are produced by Comptonevents by analyzing the detected information, calculating the angle ofCompton scattering of said detected Compton event as a function ofenergies E1 and E2 produced by a detected Compton event and anestimation of the energy E₀ of said photons emanating from the subjectand producing an image as a function of the location informationproduced by each detected Compton event and the corresponding calculatedCompton scattering angle to produce said image.

The foregoing and other objects and advantages of the invention willappear from the following description. In the description, reference ismade to the accompanying drawings which form a part hereof, and in whichthere is shown by way of illustration a preferred embodiment of theinvention. Such embodiment does not necessarily represent the full scopeof the invention, however, and reference is made therefore to the claimsherein for interpreting the scope of the invention.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

FIG. 1 is a schematic representation of a Compton-scatter aperture whichemploys two detectors;

FIG. 2 is a schematic representation of a Compton-scatter aperture whichemploys three detectors;

FIG. 3 is a pictorial view of a nuclear imaging tomographic scanningsystem which employs the present invention;

FIG. 4 is a schematic representation of the camera which forms part ofthe scanning system of FIG. 3;

FIGS. 5-7 are graphic representations showing the advantages of thepresent invention when used in the scanning system of FIG. 3 withincident photon energies of 140, 80 and 364 keV, respectively; and

FIG. 8 is a block diagram illustrating hardware for implementing theinventive method.

DETAILED DESCRIPTION OF THE INVENTION

Referring now to the drawings and specifically to FIG. 1, two detectorsD1 and D2 are positioned in parallel relationship to one side of asource S. detector D1 between source S and detector D2. A gamma particleproduced by source S having an initial energy E₀ travels along the pathindicated by arrow 2 impacting detector D1 at point X₁. Upon impact,some of the particles energy E₁ is absorbed by detector D1 and theparticle is deflected exiting detector D1 along the path indicated byarrow 4 toward detector D2 with the remaining energy E₂ (i.e. E₀−E₁).The particle impacts detector D2 at point X₂. A scatter angle betweenarrows 2 and 4 is indicated by symbol φ.

For the purposes of this explanation the measurements of the photonenergies recorded in detectors D1 and D2 will be referred to by symbolsY₁ and Y₂, respectively. Assuming measurements Y₁ and Y₂ areindependent, a posterior probability density for the deposited energiesE₁ and E₂ given measurements Y₁ and Y₂ can be expressed as:

f(E₁,E₂|Y₁,Y₂)=C×f(Y₁|E₁) f(Y₂|E₂)f(E₁,E₂|E₀)  (2)

where C is a normalization constant.

Because energy is conserved upon detection, the following generalconservation of energy equation can be written:

E₀=E₁+E₂  (3)

Equations (2) and (3) can be combined to yield:

f(E₁,E₂|Y₁,Y₂)=C×f(Y₁|E₁) f(Y₂|E₀−E₁)f(E₁|E₀)f(E₀)  (4)

Note that the prior density on energies E₁ and E₂ in Equation 2 has beensplit into the product of two terms. The first term describes theconditional probability of depositing energy E₁ in detector D1 givenincident energy E₀. Its relation to the Klein-Nishina formula (seeEquation 7 below) will become clear below. It also depends to a lesserextent on the geometry of detectors D1 and D2 and the detectorfields-of-view. For example, for a point source and small detectors,only one scattering angle would be possible, and f(E₁|E₀) would reflectthis. The prior density on the incident energies f(E₀) is simply theincident energy spectrum whose area has been normalized to unity.

In an equivalent expression the likelihood function can be written interms of the incident energy E₀ and scatter angle φ in place of energiesE₁ and E₂ as:

f(φ,E₀|Y₁,Y₂)=C×f(Y₁|E₀−g (φ,E₀))f(Y₂|g(φ,E₀))f(φ|E₀)f(E₀);  (5)

where: $\begin{matrix}{{g\left( {\phi,E_{0}} \right)} = \frac{E_{0}}{1 + {\frac{E_{0}}{m_{0}c^{2}}\left( {1 - {\cos \quad \phi}} \right)}}} & (6)\end{matrix}$

and where m₀c² is the rest mass of an election (511 keV) and f(φ|E₀) isessentially the Klein-Nishina formula normalized as a probabilitydensity function: $\begin{matrix}\begin{matrix}{{f\left( {\phi E_{0}} \right)} = \quad {C \times \sin \quad \phi \times \left( \frac{1}{1 + {\alpha\left( {1 - {\cos \quad \phi}} \right.}} \right)^{2}\left( \frac{1 + {\cos^{2}\phi}}{2} \right) \times}} \\{\quad {\left( {1 + \frac{{\alpha^{2}\left( {1 - {\cos \quad \phi}} \right)}^{2}}{\left( {1 + {\cos^{2}\phi}} \right)\left( {1 + {\alpha \left( {1 - {\cos \quad \phi}} \right)}} \right.}} \right);}}\end{matrix} & (7)\end{matrix}$

in which: $\begin{matrix}{\alpha = \frac{E_{0}}{m_{0}c^{2}}} & (8)\end{matrix}$

As noted above, f(E₁|E₀) may deviate somewhat from this formula becauseof system geometry.

Given either of the posterior density functions expressed in Equations 4or 5, an “optimum” estimator for the incident energy and scatteringangle can be derived.

The maximum a posteriori estimate for energies E₁ and E₀ (equivalent asnoted above to estimating φ and E₀) is: $\begin{matrix}\begin{matrix}{\left( {{\hat{E}}_{0},{\hat{E}}_{1}} \right) = \quad {{\underset{E_{0},E_{1}}{argmax}\text{log}{f\left( {Y_{1}E_{1}} \right)}} +}} \\{\left. \quad {{\log \quad f\left( Y_{2} \right.}{E_{0} - E_{1}}} \right) +} \\{\left. \quad {{\log\left( E_{1} \right.}E_{0}} \right) + {\log \quad {f\left( E_{0} \right)}}}\end{matrix} & (9)\end{matrix}$

If f(E₀) describes a continuous spectrum (i.e. f(E₀) contains nocharacteristic energy lines), Equation 9 can be solved by setting thegradient (with respect to E₁ and E₀) equal to zero with the appropriatechoice of estimates Ê₁ and Ê₀ such that:

∇ log f(Y₁|Ê₁)+∇ log f(Y₂|Ê₀−Ê₁) +∇ log f(Ê₁|Ê₀)+∇ log f(Ê₀)=0.  (10)

MAP estimation may be important in γ-ray astronomy applications wherethe incident spectrum can be measured with high accuracy and is sharplypeaked at several locations.

The limiting variance for estimating Ê₁ and Ê₀ for any unbiasedestimator is given by the global Cramèr-Rao (CR) bound as:

K(E₀,E₁)≧[−E[∇²(log f(Y₁,Y₂|E₀, E₁)+log f(E₁|E₀)+log f(E₀))]]⁻¹  (11)

where the K is the covariance matrix for any estimator of E₁ and E₀. MAPestimators are important because they asymptotically (in the limit of alarge number of observations) achieve the bound. Often they perform verywell relative to the bound for far fewer observations.

Although MAP estimation provides perhaps the ultimate performance, itrequires full knowledge of the incident energy spectrum or at least aprobabilistic formulation if the incident spectrum has unknowncharacteristics. Usually, information of this detailed nature is notavailable; nevertheless, the incident spectrum often containsdistinctive features. In particular, it may consist of a series of nearδ-functions or lines having characteristic energies E^((k)). Denote thesequence of incident characteristic energy lines as{E^((k))}_(k = 1)^(K).

In addition to the characteristic or discrete spectrum there is usuallya source of photons having a continuum of energies. The compositeincident spectrum is the sum of these continuous and discrete spectrumcontributions.

Often, it is only the discrete spectrum that is significant andbackground arising from the continuous spectrum is considered anuisance. This is typical in nuclear medicine where a radionuclide mayemit γ-rays at a single energy. Due to Compton-scatter in an object(e.g. a patient), however, the flux incident on the camera contains bothphotons having the characteristic energy, which are desirable, andscattered photons having a continuum of energies. Registration of thisscattered radiation results in an undesirable background. For this caseone must only decide whether or not an incident photon had an energyequal to one of the E^((k)) or not. Once an E^((k)) is selected (andthere will be some error in classifying incident photons), it can beused as a hard constraint to reduce uncertainty in the scattering angleestimate.

Rather than starting off with a discussion of the spectralclassification process, a simple estimator for scattering angle φ isdeveloped under the assumption that the incident energy E₀ is known. Theestimator does not use all information inherent in the measurement; andin particular, the prior probability density function (pdf) f(E₁|E₀) isignored. Under conditions normally encountered in medical imaging, theperformance of estimators using this prior pdf would only be slightlybetter than the simpler estimator described next.

Assuming E₀ is known and ignoring the prior pdf allows a log-likelihoodfunction for the observations, given energy E₁, to be written as:

log f(Y₁,Y₂|E₀,E₁)=log f(Y_(i)|E_(i))+log f(Y₂|E₀−E₁)  (12)

For the purpose of this explanation only, assume that the individualterms are gaussian-distributed with variances σ_(D1) ² and σ_(D2) ² fordetectors D1 and D2 respectively. Equivalently:

Y_(i)˜N(E_(i),σ_(Di) ²)  (13)

Assume also that the variances for detectors D1 and D2 are independentof energy. The maximum likelihood (ML) estimator for energy E₁ reducesto a simple least-squares problem as: $\begin{matrix}{{\hat{E}}_{1} = {\frac{\sigma_{D2}^{2}Y_{1}}{\sigma_{D1}^{2} + \sigma_{D2}^{2}} + {\frac{\sigma_{D1}^{2}\left( {E_{0} - Y_{2}} \right)}{\sigma_{D1}^{2} + \sigma_{D2}^{2}}.}}} & (14)\end{matrix}$

Equation 14 represents the ML estimator for gaussian distributed datawith known variances. Typically, for any distribution, a similargeneralized least-squares estimator can be derived where the weights aredetermined from the first- and second-moments of the likelihoodfunctions. The weights derived from the variances may depend on eitherthe data or the values of E₀ and E₁ making the estimator non-linear withrespect to the data. As is clear in Equation 14 when the incident energyis E₀ assumed known, measurements from both detectors D1 and D2 can beused to reduce errors in the estimate of E₁ and hence in scatteringangle φ.

To compare the performance of an estimator using Equation 14 withperformance of a typical Compton camera, equations which representapproximations for the variances σ_(φ) ² in scattering angle φ given theuncertainties in estimates Ê₀ and Ê₁ can be derived.

An equation for estimating the scattering angle φ from estimates ofincident energy E₀ and energy E₁ absorbed by detector D₁ upon scatteringis: $\begin{matrix}{\phi = {\cos^{- 1}\left( {1 + \frac{\alpha}{E_{0}} - \frac{\alpha}{\left( {E_{0} - E_{1}} \right)}} \right)}} & (15)\end{matrix}$

This formula is valid when Doppler broadening of the Compton line widthis inconsequential. Doppler broadening is explained below. TypicalCompton cameras have used less than all information gathered to estimatethe scattering angle φ. For example, some cameras estimate thescattering angle using either the energy measurements E₁ and E₂ from thefirst and second detectors D₁ and D₂ while others use the incidentenergy E₀ and the energy measurement from one of the two detectors D₁ orD₂. While these cameras appear to be simple and therefore desirable,unfortunately the variance in the estimates generated by these camerasis relatively large.

The variance of the scattering angle estimate is approximated bylinearizing the estimator for scattering angle and calculating thevariance of the linearized estimator.

In the first case, where energies E₁ and E₂ are measured, E₀ is unknownand the prior pdf is neglected the energies deposited in detectors D₁and D₂ can be estimated individually, which results in the followingexpression for variance σ₁₀₀ ², written in terms of E₀ and E₁:$\begin{matrix}{\sigma_{\phi}^{2} \approx {\frac{\alpha^{2}}{\sin^{2}\phi}\left\lbrack {{E_{0}^{- 4}{\sigma_{D1}^{2}\left( E_{1} \right)}} + {\left\lbrack {\frac{1}{\left( {E_{0} - E_{1}} \right)^{2}} - \frac{1}{E_{0}^{2}}} \right\rbrack^{2}{\sigma_{D2}^{2}\left( {E_{0} - E_{1}} \right)}}} \right\rbrack}} & (16)\end{matrix}$

The fact that the variances can depend on the deposited energy has beendenoted by σ_(D) ²=σ_(D) ²(E).

In the second case where E₀ is known and E₁ is measured, the variance isas follows: $\begin{matrix}{\sigma_{\phi}^{2} \approx {\frac{\alpha^{2}}{\sin^{2}\phi}\left( {E_{0} - E_{1}} \right)^{- 4}{\sigma_{D1}^{2}\left( E_{1} \right)}}} & (17)\end{matrix}$

Alternatively, in the second case where E₀ is known and E₂ is measured,the variance is expressed as: $\begin{matrix}{\sigma_{\phi}^{2} \approx {\frac{\alpha^{2}}{\sin^{2}\phi}{E_{2}}^{- 4}{\sigma_{D2}^{2}\left( E_{2} \right)}}} & (18)\end{matrix}$

These expressions for angle variance σ_(φ) ² are well known in theCompton camera art.

When the present invention is used to improve the angular estimation andenergy E₀ is known in advance, the energy measurements E₁ and E₂ fromboth detectors D1 and D2 can be used to reduce angular uncertainty basedon the conservation of energy constraint (i.e. see Equation 3). To thisend, the angular variance σ_(φ) ² using the present invention isexpressed as follows: $\begin{matrix}{\sigma_{\phi}^{2} \approx {\frac{\alpha^{2}}{\sin^{2}\phi}{{\left( {E_{0} - E_{1}} \right)^{- 4}\left\lbrack \frac{{\sigma_{D1}^{2}\left( E_{1} \right)}{\sigma_{D2}^{2}\left( {E_{0} - E_{1}} \right)}}{{\sigma_{D1}^{2}\left( E_{1} \right)} + {\sigma_{D2}^{2}\left( {E_{0} - E_{1}} \right)}} \right\rbrack}.}}} & (19)\end{matrix}$

The angular uncertainty always decreases when using two energymeasurements since: $\begin{matrix}{\frac{{\sigma_{D1}^{2}\left( E_{1} \right)}{\sigma_{D2}^{2}\left( {E_{0} - E_{1}} \right)}}{{\sigma_{D1}^{2}\left( E_{1} \right)} + {\sigma_{D2}^{2}\left( {E_{0} - E_{1}} \right)}} < {{\sigma_{D1}^{2}\left( E_{1} \right)}\quad {and}\quad {\sigma_{D1}^{2}\left( {E_{0} - E_{1}} \right)}}} & (20)\end{matrix}$

As proof, by setting a value β equal to σ_(D1) ²(E1)/σ_(D2) ²(E₀−E₁) (orvice-versa for detector D2) in Equation 20 and simplifying, we get:$\begin{matrix}{{\frac{1}{1 + \beta}\sigma_{D1}^{2}} < \sigma_{D1}^{2}} & (21)\end{matrix}$

which holds for all β>0 (the only possibility) because (1β)⁻¹<1.

It is apparent that the variance improves the most over the best singlemeasurement when σ_(D1) ²=σ_(D2) ². At this point the variance in thescattering angle is reduced by a factor of two, which results inresolution gains of 2 and 2{square root over (2)} in area and volumeresolution, respectively. Note that if the energy resolution of onedetector is much worse than that of the other as is the case in currentCompton cameras, there will be virtually no improvement in angularuncertainty. However, it is a teaching of the present invention that byusing modern solid-state detectors, it is reasonable to constructCompton-scatter apertures using two detectors having excellent energyresolution. One can thus achieve variance reductions closer to themaximum achievable factor of two for lower energy photons (<150-200keV).

Referring to FIG. 3, there is shown generally at 10, a nuclear imagingtomographic scanning system which includes a tomographic scanner 11, anda patient support table 12. The scanner 11 comprises an annular gantry13 supported in a vertical position as shown by a pedestal 14 and havinga camera 16 supported from the gantry 13 in cantilevered fashion by anarm assembly 17 and balanced by a counterweight 18 on the other end ofthe arm assembly 17. The arm assembly 17 is so connected to the gantry13 as to allow the entire arm assembly 17 to be rotated within thegantry 13 by a motor-drive system (not shown), to thereby rotate thecamera 16 in a circular path to a variety of view angles around thepatient 19 supported on the table 12. The movement of the camera 16allows the collection of multiple views which can be used to reconstructa tomographic image of the patient in the area of concern. The structureand operational movement of the scanner 11 is of a conventional nature.

Referring to FIGS. 3 and 4, as is well known in the art, the variousisotopes used in nuclear medicine are preferentially absorbed by thepatient 19 to emit gamma ray photons in a pattern which permitsvisualizing the configuration of body tissue and blood vessels. Thecamera 16 used to detect and locate the source of such emissionsincludes a first detector array D1 and a second detector D2 which abutsthe outer perimeter of the first detector D1 and extends perpendiculartherefrom.

A gamma ray emanating from a source location indicated at 22 and havingenergy E₀ passes into detector D1 and produces a Compton event at alocation X₁ in the first detector D1. The gamma ray gives up energy E₁during the Compton event and a gamma ray is emitted with energy E₂ fromlocation X₁ which is absorbed at a location X₂ in the second detectorD2. Detectors D1 and D2 produce signals indicative of sensed events, thesignals provided to acquisition circuit 28 via data buses 24 and 25,respectively.

To generate required data (i.e. time of impact, energy absorbed,location of impact), each detector D1, D2 includes means for determining(1) the arrival time of each interaction T₁ and T₂, respectively, (2)the location X₁, X₂, respectively, of each interaction within thedetector, and (3) the energy E₁, E₂, respectively, deposited during eachinteraction.

Detectors D1 and D2 preferably have similar energy resolutions. In apreferred configuration, scattering detector D1 is constructed of amaterial of low atomic number, which has low Doppler-broadening and ahigh probability that gamma-ray interactions will be Compton-scatteringrather than photoelectric absorption or coherent scattering. Seconddetector D2 is constructed of a material having higher atomic number toincrease the probability of photoelectrically absorbing the scatteredradiation. In a typical situation, scattering detector D1 is constructedusing a position-sensitive silicon detector while absorbing detector D2employs germanium or cadmium zinc telluride.

Referring to FIGS. 4 and 8, circuit 28 comprises a plurality ofcomponents including an event processor 50, an energy processor 52, ascattering angle estimator 54, a first direct image reconstructor 56 anda second image reconstructor 58 which explicitly accounts for scatterangles. Each detector D1, D2 provides an arrival time estimate T₁, orT₂, a position X₁ or X₂, and an energy measurement E₁ or E₂ for eachphoton interaction to processor 50. Processor 50 first attempts todetermine which set of detector interactions are associated with aparticular gamma-ray incident on the Compton camera. For detectorsystems exhibiting good time resolution, this association is largelyperformed as it is currently done in positron emission tomography: Iftwo (or more) interaction times fall within a short time-interval, theyare assumed to have originated from the same gamma-ray photon incidenton the Compton-camera.

The second task performed by processor 50 is to select the order ofinteractions. For example, Compton cameras are usually designed in sucha way as to enhance the probability of scattering from detector D1first. However, nothing prevents scatter from occurring in detector D2first with subsequent interaction in detector D1 (or multipleinteractions in D1 followed by D2). Processor 50 uses as its basis thelikelihood function described by Equation 32. Specifically, processor 50maximizes the likelihood term inside the integral by using estimates ofarrival times, energies and interaction positions by using time,location and energy estimates to determine order of interaction, a moreaccurate result occurs. To this end, an ordering means or calculator 69is provided in processor 50. In the alternative a suitable surrogatefunction could be used. For example, if the time resolution of eachdetector D1 and D2 is very good, the time measurements from eachdetector can be used to estimate the actual interaction times with veryhigh fidelity. In this case, processor 50 really only needs to considerthe arrival time measurements. Once the arrival times have beenestimated the interaction order follows immediately.

Processor 50 formats the events into a time-ordered list, each elementof the list corresponding to one (estimated) event incident on theCompton camera. Each list element further consists of a time-orderedsequence of detector positions and energies. For example, if the n-thgamma-ray incident on the camera interacted in detector D1 at positionX₁ with measured energy Y₁ and the scattered photon interacted indetector D2 at X₂ with energy Y₂, the n-th element of the list wouldconsist of two ordered pairs [(Y₁,X₁), (Y₂,X₂)].

Each element of the list consists of N ordered pairs where N is thenumber of detected interactions associated with a single gamma-rayincident on the Compton-camera.

The resulting event list is first processed by energy processor 52 inorder to estimate E₀. Typically, if both detectors have good energyresolution, the incident energy is estimated by first summing allmeasured energy corresponding to a single incident gamma-ray (or asingle list element) according to the following for the case of twointeractions:

E_(test)=Y₁+Y₂  (22)

This estimate of E₀ is then compared with a relatively narrow energywindow. If the event falls within the window then Ê₀ is assigned asingle energy corresponding to the characteristic energy of thegamma-ray emission of a radionuclide of interest.

As an example, consider 99m-Tc, which has a characteristic emission ofapproximately 140 keV. A narrow energy window on the Compton cameramight run from a low of 139 keV to a high of 141 keV. If E_(test) fallswithin the window then Ê₀ is assigned the value of 140 keV regardless ofthe particular value of E_(test).

Multiple energy windows can be accommodated and require multipleapplications of the method described above. If the event falls within adesired energy window, the event is augmented with the initial energyestimate Ê₀ and is passed on for further processing. If it does not fallin a desired energy range, it is removed from the list.

The list output from processor 52 may be processed in two ways. First,the scattering angles φ can be directly determined via estimator 54 fromthe measurements as described below and then passed on for imagereconstruction using reconstructor 58. Second, the list can be passeddirectly to reconstructor 56 which applies a likelihood model for thecamera system. While the second method retains more information and intheory is capable of providing better reconstructed images, the firstmethod is usually sufficient. The invention is meant to cover bothmethods of reconstructing.

Where the scattering angle φ is estimated from two interactions andDoppler broadening is negligible (assuming the interaction order isD1→D2), the following formula is used:

Ê₁=(σ₁ ²+σ₂ ²)⁻¹(σ₂ ²Y₁+σ₁ ²(Ê₀−Y₂)); and  (23)

$\begin{matrix}{\hat{\phi} = {{\cos^{- 1}\left\lbrack {\frac{\alpha}{{\hat{E}}_{0}} - \frac{\alpha}{{\hat{E}}_{0} - {\hat{E}}_{1}} + 1} \right\rbrack}.}} & (24)\end{matrix}$

Similar formulae hold for the opposite interaction sequence and for morethan two interactions. The estimated scattering angles as well as thelocation of the cone vertex (the position X₁ of the first interaction inthe Compton-camera) is added to an output list and is passed on forreconstruction to reconstructer 58.

The present invention can also be extended to cameras based onmultiple-scattering. A two-scatter aperture system is shown in FIG. 2.Assume that the geometry of the three detectors is fixed such that θ isknown. The MAP estimator can be developed in a manner analogous to thesingle-scatter case. To derive a similar estimator to that describedabove, start with the log-likelihood for observing detector outputs ofY₁, Y₂, and Y₃ given deposited energies E₁, E₂, and E₃:

 log f(Y₁,Y₂,Y₃|E₁,E₂,E₃)=log f(Y₁|E₁)+log f(Y₂|E₂)+log f(Y₃|E₃)  (25)

Assuming that E₀ is known in addition to θ provides two constraintequations:

E₃=E₀−(E₁+E₂)  (26)

and $\begin{matrix}{E_{2} = {{g_{\theta}\left( {E_{0} - E_{1}} \right)}\overset{def}{=}{\left( {E_{0} - E_{1}} \right)\left\lbrack {1 - \frac{1}{1 + {\frac{E_{0} - E_{1}}{\alpha}\left( {1 - {\cos \quad \theta}} \right)}}} \right\rbrack}}} & (27)\end{matrix}$

The corresponding constrained likelihood function is:

log f(Y₁,Y₂,Y₃|/E₁)=log f(Y₁|E₁)+log f(Y₂|g_(θ)(E₀−E₁))+logf(Y₃|E₀−(E₁+g_(θ)(E₀−E₁)))  (28)

Again, assuming the measurements are conditionally gaussian-distributed,the ML estimate of E₁ (hence the scattering angle φ) satisfies$\begin{matrix}\begin{matrix}{\frac{{\partial\log}\quad f}{\partial E_{1}} = \quad {\frac{\left( {Y_{1} - E_{1}} \right)}{\sigma_{D1}^{2}} + {\frac{\partial g_{\theta}}{\partial E_{1}}\frac{\left( {Y_{2} - {g_{\theta}\left( {E_{0} - E_{1}} \right)}} \right)}{\sigma_{D2}^{2}}} -}} \\{\quad {{\left( {1 + \frac{\partial g_{\theta}}{\partial E_{1}}} \right)\frac{\left( {Y_{3} - E_{0} + E_{1} + {g_{\theta}\left( {E_{0} - E_{1}} \right)}} \right)}{\sigma_{D3}^{2}}} = 0}}\end{matrix} & (29)\end{matrix}$

For the case where Doppler broadening is significant, there is no longera one-to-one correspondence between the scattering angle and theincident energy and energy-loss in the first interaction. In thissituation the scattering angle can be estimated by maximizing theappropriate likelihood function using Ê₀ and Ê₁ as estimated above andsolving the following equation: $\begin{matrix}{\hat{\phi} = {\arg \quad {\max\limits_{\phi}\quad {{f\left( {{\hat{E}}_{0},{{\hat{E}}_{1}\phi}} \right)}.}}}} & (30)\end{matrix}$

The probability density function in Equation 30 can be calculated usingmethods well-known in the art.

Methods for reconstructing images from a Compton camera are well knownwith the better methods based upon maximizing a penalized likelihoodobjective (sometimes referred to as Bayesian methods). All methodseither explicitly or implicitly require a model for the recordedmeasurements as a function of the underlying distribution of radiotracerin the object. With the appropriate choice of model, Compton camerareconstructions can be performed from as little as a scattering angleand cone-vertex for each event. At the other end of the spectrum, thereconstruction could process the outputs from the detectors directly tocreate an image based on the likelihood function of Equation 32 below.

Most of the estimates of scattering angle above assume that thegamma-ray Compton-scatters from a free-electron at rest. Knowing theenergy deposited in the first interaction is equivalent to knowing thescattering angle φ. As the incident gamma-ray energies decrease, thisassumption becomes poorer. Electrons in the scattering detector, whichare bound to individual atoms, possess considerable momentum. Thescattering of gamma-rays from these electrons, whose momentum is unknownpriori, is no longer described by Equation 15. Rather it is moreaccurate to assume there is a probabilistic relationship between angle φand energy E₁. Denote the probability density function (pdf) describingthis relationship as f(φ, E₁|E₀). The pdf expresses the additionaluncertainty in scattering angle due to the unknown electron motion andenergy. This uncertainty results in what has previously been referred toas Doppler-broadening of the Compton line-width.

If the order in which the scattering occurs (for example, D1→D2→D3,etc.) is known, the likelihood function for observing the measurementsgiven the unknown parameters of the incoming gamma-ray (direction andenergy) can be written as: $\begin{matrix}{{{f\left( {\phi,E_{0},{X_{1}\left\{ {Y_{1},\ldots \quad,Y_{N}} \right\}}} \right)} = {C \times {\int_{ɛ_{1}}{\cdots {\int_{ɛ_{N - 1}}{\prod\limits_{i = 1}^{N}\quad {{f\left( {Y_{i}E_{i}} \right)}{f\left( \left\{ {E_{1},\ldots \quad,E_{N - 1},{E_{0} - {\sum\limits_{j = 1}^{N - 1}\quad E_{j}}}} \right\}  \right.}E_{0}}}}}}}},{\left. \phi \right){f\left( {E_{0},\phi} \right)}{E_{1}}\quad \ldots \quad {E_{N - 1}}}} & (31)\end{matrix}$

where X₁ now denotes the location of the first interaction or the vertexof the cone of uncertainty from this relation. Equation 2 can beobtained as a special case when the relation between φ, E₁, and E₀ isdeterministic. Note that we have directly specified the likelihood forany number of scatters. In this case the scattering angle can bedetermined by maximizing Equation 31 given the measured data Y_(i) oralternatively by simply maximizing the equation:$f\left( {{\left\{ {E_{1},\ldots \quad,E_{N - 1},{E_{0} - {\sum\limits_{j = 1}^{N - 1}\quad E_{j}}}} \right\} E_{0}},\phi} \right)$

given estimates of E₀ and E₁, . . . , E_(N−1).

Generally, in the presence of this Doppler-broadening, the sequence inwhich the detectors are struck may not be as obvious as when thisdegradation is not present. The sequence can also be modeled in thelikelihood function as follows. $\begin{matrix}\begin{matrix}{{f\left( {\phi,E_{0},{X_{1}\left\{ Y_{1} \right\}},\left\{ t_{i} \right\},\left\{ X_{i} \right\}} \right)} = {C \times {\int{\cdots {\int{\left\lbrack {\prod\limits_{i = 1}^{N}\quad {{f\left( {Y_{i}E_{i}} \right)}{f\left( {t_{i}\tau_{i}} \right)}}} \right\rbrack {f\left( {\left\{ E_{i} \right\},\left\{ \tau_{i} \right\},{\left\{ X_{i} \right\} E_{0}},\phi} \right)}{f\left( {E_{0},\phi} \right)}{E_{i}}{\tau_{i}}}}}}}} & \quad\end{matrix} & (32)\end{matrix}$

where Y_(i) and t_(i) are the measured energies and arrival times,respectively. The true arrival times, energies, and positions arerepresented by τ_(i), E_(i), and X_(i). The terms inside the squarebrackets model the degradations due to the measurement system (e.g.,detector noise). Although, we have only included energy and time, wecould just as easily include measurement noise in the interactionposition.

The second conditional probability density function is determined by thegeometry of the Compton camera and by physical laws (includingdegradations such as Doppler broadening). The final term describes anyprior knowledge regarding the spectrum of incident energies andscattering angles.

Note that in this formulation the positions of the interactions as wellas the measured arrival times can be used to determine the appropriatesequence. As an example, the measured arrival times, positions andenergies can be used to estimate the sequence of true interaction times{τ₁, . . . , τ_(N)} from which it is trivial to estimate the interactionorder. The Event Processor described above can incorporate models basedon the general likelihood formulation (6) in order to improve theclassification of desired events (i.e., correct interaction order).

FIGS. 5-7 show plots of the three cases as a function of the scatteringangle for incident photon energies of 140, 80, and 364 keV. In each ofFIGS. 5 through 7 four curves are illustrated. A first curve E_(ou1)represents the case where E₀ is initially unknown which results when theenergy resolutions of detectors D1 of D2 are 1 keV FWHM and 10%,respectively. A second curve E_(ou2) represents the case where E₀ isinitially unknown when two high resolution detectors each having a 1 keVresolution, are used to determine scatter angle φ. A third curve E_(ok1)represents the case where E₀ is initially known and an energymeasurement from only one (i.e. D1 or D2) of detectors D1 or D2, wherethe detector resolution is 1 keV, is used to determine scatter angle φ.Fourth curve E_(ok2) is related to the present invention and representsthe case where E₀ is initially known, energy resolution of bothdetectors D1 and D2 is high and equal and energy measurements from bothdetectors D1 and D2 are used to determine scatter angle.

Clearly, where E₀ is unknown and an energy measurement from only asingle detector is used to determine scatter angle (i.e. curve E_(ou1))only limited performance results. When E₀ is known and only a singledetector measurement is used (i.e. curve E_(ok1)), improved performanceresults. Assuming that the energy resolutions of both detectors areequal, and combining the measurements using Equation 14 gives theultimate performance (i.e. curve E_(ok2)). When two high resolutiondetectors are used and the E₀ is not known in advance, performance (i.e.curve E_(ou2)) lies between that of the best single measurement and theoptimum method. Apparently, this method can also take advantage of theimproved resolution of both detectors. Nevertheless, at lower energies(80 and 140 keV) performance of the inventive method is significantlybetter over the angular range in which much of the scattering occurs(e.g. 0 to 90°).

To separate regions of the spectrum that contain the characteristicenergies from continuum radiation, a narrow energy window is placed onan appropriately corrected sum of the detector outputs. The better theenergy resolution of the sum of the two detectors, the better thisseparation can be accomplished.

Nothing is perfect, and inevitably some continuum radiation willcontaminate the measurements of the characteristic energies. Since thisradiation is by and large smoothly varying with respect to energy,however, its contribution to the characteristic energy windows can beestimated by placing windows that are somewhat wider to either side ofthese narrow windows. The contribution of the continuum radiation to thenarrow window is estimated, and then either explicitly or implicitlyremoved. Wide windows decrease counting noise in the measurements at theexpense of potential inaccuracy.

The corrected value from the k-th window is an estimate of the intensityof characteristic radiation having energy E^((k)). We assume thatmultiple characteristic energies are separated by intervalssignificantly larger than the energy resolution of the system so thatcross-talk will be negligible.

At higher energies a fraction of the interacting photons may not becompletely absorbed resulting in an energy of less than E₀ depositedwithin the camera. There are two alternatives for processing: either onecan not register incomplete absorptions, selecting only those eventsthat fall into the narrow window; or one can attempt to identify theselower energy events, which will now have a continuum of energies. In thefirst case, counting efficiency is reduced resulting in greater imagenoise. In the second, it may be difficult to separate desirable eventsfrom those that are not. Although an optimum separation can beaccomplished in principle using a complete model for the likelihoodfunction of the camera (Equation 32), performance may be highlysusceptible to model inaccuracies.

It should be understood that the methods described above are onlyexemplary and do not limit the scope of this invention and variousmodifications could be made by those skilled in the art that would fallunder the scope of the invention. To apprise the public of the scope ofthe this invention, we make the following claims:

What is claimed is:
 1. A camera for creating an image of the radiationdensity of a source of photons located in a subject, the cameracomprising: a first detector (D1) for detecting the time, location (X1)and energy (E1) of collisions between said first detector (D1) andphotons emanating from the subject; a second detector (D2) for detectingthe time, location (X2) and energy (E2) of collisions between saidsecond detector (D2) and photons emanating from the first detector (D1);Compton event detector means connected to the first and second detectors(D1 and D2) and being responsive to the information produced thereby toidentify sets of collisions in the respective first and second detectorsthat are produced by Compton events; and image reconstruction meansresponsive to energies E1 and E2 produced by a detected Compton eventand an estimation Ê₀ of a photon incident energy for constructing animage; angle determining means responsive to the energies E1 and E2produced by a detected Compton event and an estimation of the energy E₀of said photons emanating from the subject for calculating the angle ofCompton scattering of said detected Compton event.
 2. The camera asrecited in claim 1 which includes a third detector (D3) for detectingthe time, location and energy (E3) of collisions between said thirddetector (D3) and photons emanating from the second detector (D2); andthe angle determining means is also responsive to the energy (E3) tocalculate the angle of Compton scattering.
 3. The camera as recited inclaim 1 in which the angle determining means calculates the angle φ ofCompton scattering according to the expression:$\left( {\hat{\phi},{\hat{E}}_{0},{\hat{X}}_{1}} \right) = {\arg \quad {\max\limits_{({\phi,E_{0},X_{1}})}{f\left( {\phi,E_{0},{X_{1}\left\{ Y_{i} \right\}},\left\{ t_{i} \right\},\left\{ X_{i} \right\}} \right)}}}$

where X₁ is the point of impact on the first camera andf(φ,E₀,X_(i)|{Y_(i)},{t_(i)},{X_(i)}) is given by equation:${{C \times {\int\quad {\ldots \quad {\int\left\lbrack {\prod\limits_{i = 1}^{N}{{f\left( {Y_{i}E_{i}} \right)}{f\left( {t_{i}\tau_{i}} \right)}}} \right\rbrack}}}}}{f\left( {\left\{ E_{i} \right\},\left\{ \tau_{i} \right\},{\left\{ X_{i} \right\} E_{0}},\phi} \right)}{f\left( {E_{0},\phi} \right)}{E_{i}}{\tau_{i}}$

where C is a constant.
 4. The camera of claim 1 wherein the angledetermining means chooses an incident energy estimate E₀ as Y₁+Y₂ whereY₁ and Y₂ are measured energies corresponding to E₁ and E₂ and thedetermining means calculates the angle of the Compton scattering byfirst estimating energy E₁ using energies Y₁ and Y₂ and estimate E₀ andthen using the E₁ estimate to determine the scattering angle φ.
 5. Thecamera of claim 4 wherein the determining means estimates energy E₁ bysolving the following equation: Ê₁=(σ_(D1) ²+σ_(D2) ²)⁻¹[σ_(D2)²Y₁+σ_(D1) ²(Ê₀−Y₂)] where σ_(D1) ² and σ_(D2) ² are variances in theenergy uncertainties for the first and second detectors, respectively.6. The camera of claim 5 wherein Doppler broadening is negligible andwherein the determining means estimates angle φ by solving the followingequation:$\hat{\phi} + {\cos^{- 1}\left\lbrack {\frac{\alpha}{{\hat{E}}_{0}} - \frac{\alpha}{{\hat{E}}_{0} - {\hat{E}}_{1}} + 1} \right\rbrack}$where: $\alpha = \frac{E_{0}}{m_{0}c^{2}}$

and where m₀c² is the rest mass of an election.
 7. The camera of claim 5wherein Doppler broadening is significant and the determining meanscalculates the angle of the Compton scattering by solving the followingequation:$\hat{\phi} = {\arg \quad {\max\limits_{\phi}{{f\left( {{\hat{E}}_{0},{{\hat{E}}_{1}\phi}} \right)}.}}}$


8. A camera for creating an image of the radiation density of a sourceof photons located in a subject, the camera comprising: a first detector(D1) for detecting the time, location and energy (E1) of collisionsbetween said first detector (D1) and photons emanating from the subject;a second detector (D2) for detecting the time, location and energy (E2)of collisions between said second detector (D2) and photons emanatingfrom the first detector (D1); Compton event detector means connected tothe first and second detectors (D1 and D2) and being responsive to theinformation produced thereby to identify pairs of collisions in therespective first and second detectors that are produced by Comptonevents; image reconstruction means responsive to the energies E1 and E2produced by a detected Compton event and an estimation of the energy E₀of said photons emanating from the subject for forming an image;whereby, by being responsive to each of energies E1, E2 and E₀, imageaccuracy is increased as photon scattering angles related to energiesE1, E2 and E₀ and implicit in the image are relatively more accurate. 9.The camera of claim 8 wherein the image reconstruction means furtherincludes an angle determining means responsive to the energies E1, E2and E₀ for calculating the angle of Compton scattering of said detectedCompton event, the reconstruction means responsive to the scatteringangle for each event when producing the image.
 10. The camera of claim 8wherein the event detector also includes an ordering means fordetermining the order of detector collisions for a Compton event whereinthe ordering means determines the order of events as a function of thetime, location and energy of each collision.
 11. A method for use with acamera for creating an image of the radiation density of a source ofphotons located in a subject, the camera including first and secondphoton detectors, the method comprising the steps of: detecting thetime, location and energy (E1) of collisions between said first detector(D1) and photons emanating from the subject; detecting the time,location and energy (E2) of collisions between said second detector (D2)and photons emanating from the first detector (D1); identify pairs ofcollisions in the respective first and second detectors that areproduced by Compton events by analyzing the detected information;calculating the angle of Compton scattering of said detected Comptonevent as a function of energies E1 and E2 produced by a detected Comptonevent and an estimation of the energy E₀ of said photons emanating fromthe subject; and producing an image as a function of the locationinformation produced by each detected Compton event and thecorresponding calculated Compton scattering angle to produce said image.